#comparison

#library packages
library(ggplot2)
library(ggpubr)
library(data.table)

library(stringr)
library(optparse)

option_list <- list(
  make_option("--h", default = "", type = "character", help = "input hbulk expr file"),
  make_option("--m", default = "", type = "character", help = "input mbulk expr file"),
  make_option("--hm", default = "", type = "character", help = "input hsa meta file"),
  make_option("--mm", default = "", type = "character", help = "input mmu meta file"),
  make_option("--g", default = "", type = "character", help = "gene name"),
  make_option("--s", default = "", type = "character", help = "species"),
  make_option("--hap", default = "", type = "character", help = "hsa a group phenotypes"),
  make_option("--hbp", default = "", type = "character", help = "hsa b group phenotypes"),
  make_option("--had", default = "", type = "character", help = "hsa a group datasets"),
  make_option("--hbd", default = "", type = "character", help = "hsa b group datasets"),
  make_option("--t", default = "", type = "character", help = "test method"),
  make_option("--l", default = "", type = "character", help = "log"),
  make_option("--map", default = "", type = "character", help = "Mma a group phenotypes"),
  make_option("--mbp", default = "", type = "character", help = "Mma b group phenotypes"),
  make_option("--mad", default = "", type = "character", help = "Mma a group datasets"),
  make_option("--mbd", default = "", type = "character", help = "MMa b group datasets"),
  make_option("--mam", default = "", type = "character", help = "Mma a group Phenotype models"),
  make_option("--mbm", default = "", type = "character", help = "Mma b group Phenotype models"),
  make_option("--o", default = "", type = "character", help = "output png file")
)
opt <- parse_args(OptionParser(option_list = option_list))

#import source data #修改

hbulk<-fread(opt$h, data.table = F, sep = '\t')
rownames(hbulk)<-hbulk$gene
gene <- hbulk$gene
hbulk<-hbulk[,-1]
hmeta<-read.table(opt$hm, header = T, stringsAsFactors = F, sep = '\t')
mbulk<-fread(opt$m, data.table = F, sep = '\t')
rownames(mbulk)<-mbulk$gene
mbulk<-mbulk[,-1]
mmeta<-read.table(opt$mm, header = T, stringsAsFactors = F, sep = '\t')

#dim parameters
species<-opt$s #choose one from c('Hsa','Mmu')
test<-opt$t #choose one from c('wilcox.test','t.test')
log.scale<-opt$l #Yes or No
if (species == 'Hsa'){
  meta<-hmeta
  expr<-hbulk
  haps <- opt$hap %>%
    str_split(";") %>%
    unlist()
  hbps <- opt$hbp %>%
    str_split(";") %>%
    unlist()
  hads <- opt$had %>%
    str_split(";") %>%
    unlist()
  hbds <- opt$hbd %>%
    str_split(";") %>%
    unlist()
  group.a<-list(phenotype=haps, dataset=hads) #网页多选
  group.b<-list(phenotype=hbps, dataset=hbds) #网页多选
  meta.a<-meta[meta$Phenotype %in% group.a$phenotype & meta$Dataset %in% group.a$dataset,]
  meta.b<-meta[meta$Phenotype %in% group.b$phenotype & meta$Dataset %in% group.b$dataset,]
} else { #修改
  meta<-mmeta
  expr<-mbulk
  maps <- opt$map %>%
    str_split(";") %>%
    unlist()
  mbps <- opt$mbp %>%
    str_split(";") %>%
    unlist()
  mads <- opt$mad %>%
    str_split(";") %>%
    unlist()
  mbds <- opt$mbd %>%
    str_split(";") %>%
    unlist()
  mams <- opt$mam %>%
    str_split(",") %>%
    unlist()
  mbms <- opt$mbm %>%
    str_split(",") %>%
    unlist()
  group.a<-list(phenotype=maps, dataset=mads, model=mams)
  group.b<-list(phenotype=mbps, dataset=mbds, model=mbms)
  meta.a<-meta[meta$Phenotype %in% group.a$phenotype & meta$Dataset %in% group.a$dataset & meta$Model_show %in% group.a$model,]
  meta.b<-meta[meta$Phenotype %in% group.b$phenotype & meta$Dataset %in% group.b$dataset & meta$Model_show %in% group.b$model,]
}

#extract data

#check gene symbol
if (gene %in% rownames(expr)){
  expr.a<-expr[gene,meta.a$Run]
  expr.b<-expr[gene,meta.b$Run]
  meta.a<-mutate(meta.a,group='Group A',expr= as.numeric(expr.a),log2expr = log2(as.numeric(expr.a)+1))
  meta.b<-mutate(meta.b,group='Group B',expr= as.numeric(expr[gene,meta.b$Run]),log2expr = log2(as.numeric(expr.b)+1))
  meta.filter<-rbind(meta.a,meta.b)
} else {
  print("The gene does not exist in the samples selected. Please check if the species is chosen correctly.")
}


if (log.scale=='Yes'){
  png(file=opt$o,width=2000,height=2000,res=300)
  p<-ggplot(meta.filter,aes(x=group,y=log2expr))+geom_boxplot(linetype='dashed',outlier.size = 0.8)+
    stat_boxplot(aes(ymin = ..lower..,ymax=..upper..,fill=group),outlier.colour=NA)+
    scale_fill_manual(values = c("#FB8072","#80B1D3"))+
    stat_boxplot(geom = 'errorbar',aes(ymin = ..ymax..),width=0.3)+
    stat_boxplot(geom = 'errorbar',aes(ymax = ..ymin..),width=0.3)+
    labs(y='Log2(TPM+1)',title = gene)+ #修改
    theme_bw()+theme(legend.position = 'none',axis.title.y = element_text(face = 'bold'),axis.title.x = element_blank(),
                     plot.title = element_text(hjust = 0.5,face = 'bold'))+
    stat_compare_means(method = test,comparisons = list(c('Group A','Group B')),label = 'p.signif',bracket.size = 0.6)
  print(p)
  dev.off()
} else {
  png(file=opt$o,width=2000,height=2000,res=300)
  p<-ggplot(meta.filter,aes(x=group,y=expr))+geom_boxplot(linetype='dashed',outlier.size = 0.8)+
    stat_boxplot(aes(ymin = ..lower..,ymax=..upper..,fill=group),outlier.colour=NA)+
    scale_fill_manual(values = c("#FB8072","#80B1D3"))+
    stat_boxplot(geom = 'errorbar',aes(ymin = ..ymax..),width=0.3)+
    stat_boxplot(geom = 'errorbar',aes(ymax = ..ymin..),width=0.3)+
    labs(y='TPM',title = gene)+
    #scale_y_continuous(labels = function(x)format(x,scientific = F,digits = 0))+
    theme_bw()+theme(legend.position = 'none',axis.title.y = element_text(face = 'bold'),axis.title.x = element_blank(),
                     plot.title = element_text(hjust = 0.5,face = 'bold'))+
    stat_compare_means(method = test,comparisons = list(c('Group A','Group B')),label = 'p.signif',bracket.size = 0.6)
  print(p)
  dev.off()
}

#export csv
write.csv(meta.filter,file = paste('gene.csv',sep = ''),row.names = F)
